#include "AlgorithmSetup.h"

#include <stdio.h>

#include "math.h"
#include "string.h"

//static const double PI = atan(1.0) * 4.0;
static const double PI = 3.141592653589;

int GetFourierWindowSize(int SpectrumSize)
{
	return SpectrumSize * 2;
}

void GenerateFourierWindow(double MaxDetectFrequency, int SpectrumSize, int SampleCount, double *TimeVector, double *WindowRe, double *WindowIm)
{
	double DeltaF = MaxDetectFrequency / SpectrumSize * 2.0;
	int i, j;
	int WindowSize = GetFourierWindowSize(SpectrumSize);

	// all Re[k] and Im[k] need to be zero.
	memset(WindowRe, 0, WindowSize * sizeof(double));
	memset(WindowIm, 0, WindowSize * sizeof(double));

	// This is essentially a NDFT for a signal with a frequency of ZERO and a phase of ZERO,
	// i. e. the real part of the signal is always 1 and the imaginay part always 0.
	// the "normal" ndft formulae still apply:
	// sum_re = sum_re + (SamplesRe[i] * cos(arg) - SamplesIm[i] * sin(arg)) / SampleCount;
	// sum_im = sum_im + (SamplesRe[i] * sin(arg) + SamplesIm[i] * cos(arg)) / SampleCount;
	// however with  SamplesRe[i] = 1  and  SamplesIm[i] = 0 , the calculation is simplified
	for (i = 0; i < WindowSize; i++)
	{
		double Frequency = DeltaF * i - MaxDetectFrequency * 2.0;
		for (j = 0; j < SampleCount; j++)
		{
			double Arg = -2.0 * PI * Frequency * TimeVector[j];
			WindowRe[i] += (double)cos(Arg);
			WindowIm[i] += (double)sin(Arg);
		}
	}
}

int GetSinCosArraySize(int SampleCount, int SpectrumSize)
{
	return SampleCount * SpectrumSize;
}

void GenerateSinCosArray(int SampleCount, double *TimeVector, double MaxDetectFrequency, int SpectrumSize, double *SinArray, double *CosArray)
{
	double DeltaF = MaxDetectFrequency / SpectrumSize * 2.0;
	int i, j;

#ifdef WIN32
#define FREQARG_TXT "freqarg_windows.txt"
#else
#define FREQARG_TXT "freqarg.txt"
#endif
	FILE *fp = fopen(FREQARG_TXT, "w");
	
	printf("DeltaF = %lf\nMaxDetectFrequency = %lf\nSpectrumSize = %d\nSampleCount = %d\n", DeltaF, MaxDetectFrequency, SpectrumSize, SampleCount);

	for (i = 0; i < SampleCount; i++)
	{
		for (j = 0; j < SpectrumSize; j++)
		{
			double Frequency = DeltaF * j - MaxDetectFrequency;
			double Arg = -2.0 * PI * Frequency * TimeVector[i];
			fprintf(fp, "%lf, %lf\n", Frequency, Arg);
			CosArray[j * SampleCount + i] = (double)cos(Arg);
			SinArray[j * SampleCount + i] = (double)sin(Arg);
		}
	}
	fclose(fp);
}






























